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1. Introduction 

In the last few years a big effort has been devoted by several authors to the problem of an 
efficient computation of one-loop corrections for multi-particle processes. This is a problem 
relevant for both LHC and ILC physics. In the case of QCD, the NLO six gluon amplitude 
has been recently obtained by two different groups [1], and, in the case of e + e~ collisions, 
complete EW calculations, involving 5-point [2] and 6-point [3] loop functions are available 
at the cross section level. The used techniques range from purely numerical methods to an- 
alytic ones, also including semi-numerical approaches. For analytical approaches, the main 
issue is reducing, using computer algebra, generic one-loop integrals into a minimal set of 
scalar integrals (and remaining pieces, the so called rational terms), mainly by tensor re- 
duction [4-7]. For multi-particle processes though this method becomes quite cumbersome 
because of the large number of terms generated and the appearance of numerical insta- 
bilities due to the zeros of Gram-determinants. On the other hand, several numerical or 
semi- numerical methods aim for a direct numerical computation of the tensor integrals [8] . 
Although purely numerical methods can in principle deal with any configuration of masses 
and also allow for a direct computation of both non-rational and rational terms, their 
applicability remains limited due to the high demand of computational resources and the 
non-existence of an efficient automation. 

In a different approach, the one-loop amplitude rather than individual integrals are 
evaluated using the unitarity cut method [9], which relies on tree amplitudes and avoids 
the computation of Feynman diagrams. In another development, the four-dimensional 
unitarity cut method has been used for the calculation of QCD amplitudes [10], using 
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twistor-based approaches [11]. Moreover, a generalization of the the unitarity cut method 
in d dimensions, has been pursued recently [12]. 

Nevertheless, in practice, only the part of the amplitude proportional to the loop 
scalar functions can be obtained straightforwardly. The remaining piece, the rational part, 
should then be reconstructed either by using a direct computation based on Feynman 
diagrams [13-15] or by using a bootstrap approach [16]. Furthermore the complexity of 
the calculation increases away from massless theories. 

In a recent paper [17], we proposed a reduction technique for arbitrary one-loop sub- 
amplitudes at the integrand level by exploiting numerically the set of kinematical equations 
for the integration momentum, that extend the quadruple, triple and double cuts used in 
the unitarity-cut method. The method requires a minimal information about the form of 
the one-loop (sub-)amplitude and therefore it is well suited for a numerical implementation. 
The method works for any set of internal and/or external masses, so that one is able to 
study the full electroweak model, without being limited to massless theories. 

In this paper, we describe our experience with the first practical non-trivial imple- 
mentation of such a method in the computation of a physical process: namely 27 — > 47, 
including massive fermion loops. For the massless case, there are a few results available 
in the literature. Analytical expressions were first presented by Mahlon [18] some time 
ago, however his results do not cover all possible helicity configurations. More recently the 
complete set of six-photon amplitudes was computed numerically by Nagy and Soper [19]. 
Very recently the same results were also obtained by Binoth et al. [20], that also provide 
compact analytical expressions. 

In section 2, we recall the basics of our method and, in particular, we show how the 
knowledge of the rational terms can be inferred, with full generality, once the coefficients 
of the loop functions have been determined. 

In section 3, we outline our solution to cure the numerical inaccuracies related to the 
appearance of zeros of Gram-determinants. We explicitly illustrate the case of 2-point 
amplitudes, that we had to implement to deal with the process at hand. 

In section 4, we present our numerical results. For massless fermion loops we compare 
with available results. Moreover, since we are not limited to massless contributions, we 
also present, for the first time, results with massive fermion loops. 

Finally, in the last section, we discuss our conclusions and future applications. 

2. The method and the computation of the rational terms 

The starting point of the method is the general expression for the integrand of a generic 
m-point one-loop (sub-)amplitude [17] 

A(q) = N{q) , A = (9+Pi) 2 -m?, (2.1) 

L> L>i ■ ■ ■ L> m -i 

where we use a bar to denote objects living in n = 4 + e dimensions, and q 2 = q 2 + q 2 
1 . In the previous equation, N(q) is the 4-dimensional part of the numerator function of 

1 q 2 is e-dimensional and (q ■ q) = 0. 
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the amplitude 2 . N(q) depends on the 4-dimensional denominators A = (q + Pi) 2 — m 2 as 
follows 

m—l m—l 

n (q) = XT d(i iii2i3) + d(q;ioiii2i3) }[ Di 

i()<n<*2<«3 i^*0,H,«2,*3 
m—l m—l 

+ XT i^ohh) + c(q;iohi 2 )] j] A 

io<ii<i2 i^»o>*l>*2 
m—l m—l 

+ [6(io«i) +6(<z;io*i)] II A 

io<ii i^io,'h 
m—l m—l 

+ [o(io)+a(g;i )] II A 

+ P(g) H A • (2.2) 

i 

Inserted back in Eq. (2.1), this expression simply states the multi-pole nature of any m- 
point one-loop amplitude, that, clearly, contains a pole for any propagator in the loop, 
thus one has terms ranging from 1 to m poles. Notice that the term with no poles, 
namely that one proportional to P(q) is polynomial and vanishes upon integration in 
dimensional regularization; therefore does not contribute to the amplitude, as it should be. 
The coefficients of the poles can be further split in two pieces. A piece that still depend 
on q (the terms d,c,b,a), that vanishes upon integration, and a piece that do not depend 
on q (the terms d, c, b, a). Such a separation is always possible, as shown in Ref. [17], and, 
with this choice, the latter set of coefficients is therefore immediately interpretable as the 
ensemble of the coefficients of all possible 4, 3, 2, 1-point one-loop functions contributing 
to the amplitude. 

Once Eq. (2.2) is established, the task of computing the one-loop amplitude is then 
reduced to the algebraical problem of determining the coefficients d, c, b, a by evaluating the 
function N(q) a sufficient number of times, at different values of q, and then inverting the 
system. That can be achieved quite efficiently by singling out particular choices of q such 
that, systematically, 4, 3, 2 or 1 among all possible denominators A vanishes. Then the 
system of equations is solved iteratively. First one determines all possible 4-point functions, 
then the 3-point functions and so on. For example, calling q^ the 2 (in general complex) 
solutions for which 

D = A = D 2 = D 3 = , (2.3) 

(there are 2 solutions because of the quadratic nature of the propagators) and since the 
functional form of d(q; 0123) is known, one directly finds the coefficient of the box diagram 
containing the above 4 denominators through the two simple equations 

N((fi) = [d(0123)+ d((^;0123)] H A(^)- (2-4) 

i^0,l,2,3 

2 If needed, the e-dimensional part of the numerator should be treated separately, as explained in [21]. 
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This algorithm also works in the case of complex denominators, namely with complex 
masses. Notice that the described procedure can be performed at the amplitude level. One 
does not need to repeat the work for all Feynman diagrams, provided their sum is known: 
we just suppose to be able to compute N(q) numerically. 

As a further point notice that, since the terms d, c, b, a still depend on q, also the 
separation among terms in Eq. (2.2) is somehow arbitrary. Terms containing a different 
numbers of denominators can be shifted from one piece to the other in Eq. (2.2), by relaxing 
the requirement that the integral over the terms containing q vanishes. This fact provides 
an handle to cure numerical instabilities occurring at exceptional phase-space points. In 
Section 3 we will show in detail such a mechanism at work for the 2-point part of the 
amplitude. 

The described procedure works without any modification in 4 dimensions. However, 
even when starting from a perfectly finite tensor integral, the tensor reduction may even- 
tually lead to integrals that need to be regularized. A typical example are the rank six 
6-point functions contributing to the scattering 27 — > 47 we want to study. Such tensors 
are finite, but tensor reduction iteratively leads to rank m m-point tensors with 1 < m < 5, 
that are ultraviolet divergent when m < 4. For this reason, we introduced, in Eq. (2.1), 
the d-dimensional denominators A, that differs by an amount q 2 from their 4-dimensional 
counterparts 



The result of this is a mismatch in the cancellation of the <i-dimensional denominators of 
Eq. (2.1) with the 4-dimensional ones of Eq. (2.2). The rational part of the amplitude 
comes from such a lack of cancellation. 

In [17] the problem of reconstructing this rational piece has been solved by looking at 
the implicit mass dependence in the coefficients d, c, b, a of the one-loop functions. Such 
a method is adequate up to 4-point functions; for higher-point functions the dependence 
becomes too complicated to be used in practice. In addition, it requires the solution 
of further systems of linear equations, slowing down the whole computation. For those 
reasons, we suggest here a different method. One starts by rewriting any denominator 
appearing in Eq. (2.1) as follows 



A = A + f ■ 



(2.5) 



A 



1 



Zi, 
A' 



with 




(2.6) 



This results in 



N(q) 



ZqZ\ ■ ■ ■ Z m -i . 



(2.7) 



A) A ■ A 



m— 1 



Then, by inserting Eq. 



(2.2) in Eq. (2.7), one obtains 
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m-1 • • • \ m-1 

+ c ^o^ ? 2j + c(g;* *i* 2 ) "Q 2. 

io<u<«2 u 1 ^ j^io,n,«2 

^ 6(i ii) + 6(9; ion) fr 7 
+ 2^ 5-5^ 11 z * 

10<11 u 1^0 ,u 

m-1 /■ \ , -/ • \ m-1 

+ a(z ) + a(g;*o) j-j ^ 

m—1 

+ P(q) H Zi . (2.8) 

i 

The rational part of the amplitude is then produced, after integrating over d n q, by the 
q 2 dependence coming from the various Zi in Eq. (2.8). It is easy to see what happens, 
for any value of m, by recalling the generic q dependence of the spurious terms. In the 
renormalizable gauge one has [17] 

P(q) = 0, 
a(q;i ) = a^(i ; \){q+p io )^ , 
6(g;i ii) = & M (ioii; + Pi )n + ^{kh; 2)(o + Pi ) ll {q + Pi )„ , 
c(q; i *i* 2 ) = c^(i *i* 2 ; l)(q + Pio)» + of*" (iohh; 2)(q + Vi )M + Pijv , 
+ ff ivp (i iiii;S){q + p iQ )^{q + Pi )„(q + Pi ) p , 
%; *'oiii 2 i 3 ) = «i M (*oii*2«3; + Pi )fi ■ ( 2 -9) 

Eq. (2.9) simply states the fact that 5(</;io) and d(q; 10*1*2*3) are at most linear in (g+Pi ), 
«o*i) at most quadratic, and c(g; 20*1*2) at most cubic. The tensors denoted by (• • • ; 1), 
(• • • ; 2) and (■ • • ; 3) stand for the respective coefficients. We will also make use of the fact 
that, due to the explicit form of the spurious terms [17] 

cf 11 '(*'o*i*2; 2) g^ = 0, 

c^ p (*o*'i* 2 ; 3) g^ = c^(i iii 2 ; 3) g w = c^(i iii 2 ; 3) g vp = and 

6^(i ii;2)^ = 0. (2.10) 

The necessary integrals that arise, after a change of variable q — > q — pi , are of the form 

T {n;2t> f rg ~ 2 £ gf£i ' ' ' %r 

D(ki) = (q + h) 2 - m 2 , h = Pi - Po (/c = 0) , (2.11) 

where we used a notation introduced in [22] and r < 3. Such integrals (from now on called 
extra-integrals) have dimensionality V = 2(1 + ^ — s) + r and give a contribution O(l) only 
when T> > 0, otherwise are of 0(e). This counting remains valid also in the presence of 
infrared and collinear divergences, as explained, for example, in Appendix B of [22] and 
in [14]. 
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We also note that, since all are a-dimensional, the dimensionality T> of the extra- 
integrals generated through Eq. (2.8) does not depend on m. We list, in the following, all 
possible contributions, collecting the computational details in Appendix A. 

Contributions proportional to G?(io«i*2*3) 

In this case r = 0. All extra-integrals are therefore scalars with V = — 4 and do not 
contribute. 

Contributions proportional to ^(^0*1^2*3; 1) 

In this case r = 1. All extra- integrals are therefore rank one tensors with V = — 3 and do 
not contribute. 

Contributions proportional to c(«oii*2) 

In this case r = with V = — 2 and no contribution 0(1) is developed. 
Contributions proportional to c^(io*i*2; 1) 

Here r = 1 and T> = — 1. Therefore, once again, there is no contribution. 
Contributions proportional to ct 1,1 ' {i§i\i2\ 2) 

Now r = 2 with P = and a finite contribution is in principle expected, generated by 
extra-integrals of the type 

4^2)) . (2 . 12) 

Nevertheless, such contribution is proportional to g^ v [22]. Therefore, due to Eq. (2.10), it 
vanishes. 

Contributions proportional to c^ vp {ioiii2; 3) 

Now r = 3 and T> = 1. The contributing extra- integrals are of the type 

&~ 2)) > ( 2 - 13 ) 

and one easily proves that the contributions 0(1) are always proportional to or <7 Mp or 
g vp . Therefore, thanks again to Eq. (2.10), they vanish. 

Contributions proportional to 6(io«i) 

Those are the first non vanishing contributions. The relevant extra-integrals have r = 
and V = 

r(n;2( S -l)) 

with 2 < s < m — 1. They have been computed, for generic values of s, in [22] (see also 
Appendix A) 

7 (n;2( S -l)) = _ i7r 2 1 + Q ,s 
s(s - 1) 
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Contributions proportional to 6 M (io«i;l) 

In this case the relevant extra-integrals are 4- vectors with V = 1 

4^ 2(S_1)) with 2 < s < m - 1 . 
A computation for generic values of s gives 



Contributions proportional to (ioii; 2) 

The relevant extra-integrals are now rank two tensor with P = 2 



i S (S_1) with 2 < s < m - 1 



They read 



/(-.2(-D) = -2Z7T 2 - 

"'^ (s + 2)(s + 

+ O(g^) + 0(e) . (2.16) 

The g^ v part is never needed because ^"(ioii; 2) g^ v = 0, according to Eq. (2.10). 
Contributions proportional to a(io) 
They involve scalar extra-integrals with V = 2 

li n]2s) , with 1< s < m - 1 . 

One computes 



1^ = -2i7r 2 - i J V fc? + I V Y( kj ■ ki) + £_lL^ y , 

V A 7 ^J=l 3=1 i+3 3=0 



+ 0{e). (2.17) 
Contributions proportional to o M (zo; 1) 

This last category involves extra-integrals with r = 1 and P = 3 

4^ 2s) , with 1< a < m - 1 . 

One obtains 

1 



2 



j(n;2s) _ 
S; " (s + 3)(s + 2)( S + l)s 
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s s s 



3=1 i^j Ij^i 



2Y J {m 2 j -k 2 j ){k j ), 

3=0 



j=0 i^j 



+ C(e). (2.18) 



To conclude, the set of the five formulas in Eqs. (2.14)-(2.18) allows one to compute the 
rational part of any one-loop m-point (sub-)amplitude, once all the coefficients of Eq. (2.2) 
have been reconstructed. 

3. Dealing with numerical instabilities 

In this section we show how to handle, in the framework of the method illustrated in the 
previous section, the simplest numerical instability appearing in any one-loop calculation, 
namely that one related to the tensor reduction of 2-point amplitudes in the limit of 
vanishing Gram-determinant 3 . This situation is simple enough to allow an easy description, 
but the outlined strategy is general and not restricted to the 2-point case. 

We start from the integrand of a generic 2-point amplitude written in the form 

in which we suppose N(q) at most quadratic in q. Our purpose is dealing with the situation 
in which k\ = (pi — po) 2 = exactly (that always occur in processes with massless external 
particles) , as well as to set up an algorithm to write down approximations around this case 
with arbitrary precision. 

According to Eq. (2.2), we can write an expansion for N(q) as follows: 

N(q) = [b(01) + b(q; 01)] + [a(0) + a(q; 0)]£>i + [a(l) + a(q; 1)]D . (3.2) 

If the Gram-determinant of the 2-point function is small, the reduction method introduced 
in [17] cannot be applied, because the solution for which Dq = D\ = 0, needed to determine 
the coefficients b and b, becomes singular 4 , in the limit of k\ — > 0, when adding the 
requirement 

J d n qb(q;0l) =0. (3.3) 

Then, we must consider two separate cases: 

k\ , but kl ^ , 

k\ -> , because h% = . (3.4) 



3 In this case the Gram-determinant is simply the square of the difference between the momenta of the 
two denominators. 

4 Such a solution goes like 1/kf. 



-8- 



The former situation may occur because of the Minkowskian metric, while the latter takes 
place at the edges of the phase-space, where some momenta become collinear. In the first 
case one can still find a solution for which Do = D\ = by relaxing the further requirement 
of Eq. (3.3). Such a solution is given in Appendix B and goes like \/{k\.v), where v is an 
arbitrary massless 4- vector, therefore is never singular in the first case of Eq. (3.4). The 
price to pay is that new non zero integrals appear of the type 5 

[(q+Po) ■ v] 



d! l q l - 



D D 1 



with j = 1, 2 and v 2 = . 



(3.5) 



What has been achieved with this new basis is then moving part of the 1-point functions 
to the 2-point sector, in such a way that combinations well behaved in the limit k\ — ► 
appear. The fact that solutions exist to the condition Dq = D\ = 0, still allows one to 
find the coefficients of such integrals (together with all the others). This solves the first 
part of the problem, namely reconstructing N(q) without knowing explicitly its analytic 
structure, but one is left with the problem of computing the new 2-point integrals. In the 
following, we present our method to determine them at any desired order in k\. Let us 
first consider the case j = 1 in Eq. (3.5). The contribution 0(1) can be easily obtained 
from the observation that 6 



/ 



d n q 



(q ■ v)(q ■ h 



0{k 2 ) 



(3.6) 



D(ko)D(ki) 

as it is evident by performing a tensor decomposition. On the other hand, by reconstructing 
denominators, one obtains 



(q-hf 



+ 



D(h) - D(ko) 



(q ■ h) + 



/ 



with 



t 2 j 2 2 



(3.7) 



(3.8) 



Eq. (3.7), inserted in Eq. (3.6) gives the desired expansion in terms of loop functions with 
less points but higher rank, in agreement with well know results [23, 24] 

(q-v) 1 f _ , . ( \ 1 



/ 



d n q- 



d n q{q-v) 



1 +^H)+0(k 2 1 ).(3.9) 



'Dik^Dih) fj ^ 'KDih) D(k )J\ f 

Expansions at arbitrary orders in k\ can be obtained in an analogous way from the two 
following equations: 

(q-hf-- 



I 



d n q 



f_\ p D^-Djko) 
2 2 



(q ■ v)(q ■ h) 2 P _ 2p 
D(ko)D( kl ) ~ U(kl h 



E 



i+j=p—l 



(q-hy 



f 



(3.10) 



5 Since v 2 = they still fulfill the third one of Eqs. (2.10), therefore, even in this case, terms 0(g l _ lv ) can 
be neglected in Eq. (2.16). 

6 From now on, we shift the integration variable: q — » q — po- The definition of the new resulting 
denominators is given in Eq. (2.11). 
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To deal with the case j = 2 in Eq. (3.5) one starts instead from the equation 

dQ D(ko)D(h) -°^>- 

This procedure breaks down when the quantity / vanishes. In this case a double expansion 
in k\ and / can still be found in terms of derivatives of one-loop scalar functions. We 
illustrate the procedure for the case j = 1 of Eq. (3.5). Our starting point is now the 
equation 

D(k )=D(k 1 )-2(q-k 1 ) + f. (3.12) 
By multiplying and dividing by D(ko) one obtains 

Applying once more Eq. (3.12) to the last integral gives 

f n (q-v)(q-h) _ f (q-vXq-h) , , 
J d q D{k ?D{ kl ) ~ J d q D{k Q fD{ kl ) [D{h) 2{q kl) + n 

Since the last integral in the previous equation is 0(kf), the final result reads 

(3.15) 

In a similar fashion, expansions at any order can be obtained. 

We now turn to the second case of Eq. (3.4), namely k^ — > 0. In this case no solution 
can be found to the double cut equation 

D(k ) = D(h) =0. (3.16) 

The reason is that now D(k\) and D(ko) are no longer independent: 

D(k ) = D(h) + f + 0(h) , (3.17) 

and clearly no q exists such that the two denominators can be simultaneously zero. Notice 
that this also implies that one cannot fit separately the coefficients of the 2-point and 1- 
point functions in Eq. (3.2). This results is a singularity l/(ki ■ v) in the system given of 
Appendix B and we should change our strategy. We than go back to Eq. (3.1) and split 
the amplitude from the beginning by multiplying it by 

u J M -fl W + %_t ! ) | (318) 
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resulting to 



A(q)=AW{q) + A<V(q) + 0(k 1 ), 



(3.19) 



with 



(3.20) 



Now the two amplitudes A^ 1 ' 2 ^ can be reconstructed separately, without any problem of 
vanishing Gram-determinant. Notice also that corrections at orders higher than 0(k\) are 
perfectly calculable by inserting again Eq. (3.18) in the term 0{k\) of Eq. (3.19). 

Once again, when / — > 0, double expansions in k\ and / can be obtained involving 
derivatives of scalar loop functions by using Eq. (3.12). For example, at the zeroth order 
in k\ and at the first one in /, one gets 



A{q) 



N{q) 



N{q) 



D(ko)D(h) Diko^Dih) 



[D(h) - 2(q ■ h) + f 



N(q) 
D(k f 
N(q) 



+ f 



+ f 



N(q) 



DikofDih 
N(q) 



[D(k 1 )-2(q-k 1 )+f]+0(k 1 ) 



D(k ) 2 J D(k 



+ 0(h) + 0(f 2 ) . 



(3.21) 



This last case exhausts all possibilities. 

The same techniques can be applied for higher-point functions. For example, in the 
case of a 3-point function, instead of k\, one introduces the 4- vector 



s^ 1 = det 



(*2 • fcl) 



(k 2 ■ k 2 



with the properties 



s ■ &2 = , s 2 oc A(ki, k 2 ) , (k\ ■ s) oc A(/ci, k 2 ) 



(3.22) 



(3.23) 



where A(/ci, k 2 ) is the Gram-determinant of the two momenta k\ and k 2 . Then, instead of 
Eq. (3.6) one has, for example, 



d n q 



(q ■ v)(q ■ s) 



L>(fco)£>(&i)£>(fc2) 



0(A(h,k 2 )) 



(3.24) 



As before, A(k\,k 2 ) can vanish either because s 2 = or / = and the two cases should 
be treated separately. 



4. Results and comparisons 

We started by checking our implementation of the rational terms. For 4-point functions up 
to rank four, we reproduced the results obtained with the alternative technique illustrated 
in [17]. Furthermore, we reproduced the rational part of the full 27 — > 2j amplitude given 
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in [25]. We also checked with an independent calculation [26] the rational terms coming 
from all of the 6-point tensors up to rank six. Finally, we computed the rational piece of 
the whole 27 — > 47 amplitude by summing up all 120 contributing Feynman diagrams and 
finding zero, as it should be [14]. 

As a first test on full amplitudes, we checked our method by reproducing the contribu- 
tion of a fermion loop to the 27 — > 27 process. This result is presented in Eqs. (A.18)-(A.20) 
of Ref. [25], for all possible helicity configurations. We are in perfect agreement with the 
analytic result, in both massless and massive cases. 

The next step was the computation of the 27 — > 47 amplitude with zero internal mass 7 , 
finding the results given in Fig. 1 and Fig. 2. It should be mentioned that our results are 
obtained algebraically, so there is no integration error involved. In Fig. 1, we reproduce 

25000 

20000 

m e 15000 
1 

03 10000 
5000 



0.5 1 1.5 2 2.5 3 

e 

Figure 1: Comparison with Fig. 5 of Ref. [19]. Helicity configurations [+ H ] and 

[H M — ] for the momenta of Eq. (4.1), represented by black dots and gray diamonds re- 
spectively, and comparison with the analytic result of Ref. [18] (continuous line). 

the results presented by Nagy and Soper [19] and very recently also by Binoth et al. [20]. 
We employ the same values of the external momenta as in Fig. 5 of Ref. [19], namely the 
following selection of final state three-momenta {p3,P4,P5,Pe}- 

p 3 = (33.5,15.9,25.0), 
p 4 = (-12.5,15.3,0.3), 
p 5 = (-10.0,-18.0,-3.3), 

Pa = (-11.0, -13.2, -22.0) . (4.1) 

After choosing the incoming photons such that they have momenta p\ and P2 along the 
z-axis, we present in the plot the amplitude obtained by rotating the final states of an- 
gle 6 about the y-axis. This is done for both helicity configurations [+ H ] and 

7 We thank Andre van Hameren for providing us with his program to compute massless one-loop scalar 
integrals. 




- 12 - 



18000 
16000 
14000 

CO 
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-2- 10000 
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e 

Figure 2: Hclicity configurations [+ H ] and [+ H 1 — ] for the momenta of Eq. (4.2), 

represented by black dots and gray diamonds respectively, and comparison with the analytic result 
of Ref. [18] (continuous line). 

[H h H — ]■ In the same plot also appears the analytic results for the configuration 

[+ H ] obtained by Mahlon [18]. In Fig. 2, we use a different set of external mo- 
menta. Starting from the following choice of {p3,P4,P5,Pe}- 

p 3 = (-10.0,-10.0,-10.0), 
p A = (12.0,-15.0,-2.0), 
p 5 = (10.0,18.0,3.0), 

p 6 = (-12.0,7.0,9.0) (4.2) 

we proceed as in the previous case. The results for the amplitudes are plotted in Fig. 2 
for the helicity configurations [+ H ] and [+ H 1 — ]. It is known that the six- 
photons amplitude vanish for the helicity configurations [+ + + + ++] and [+ + + + H — ], 
we checked this result for both choices of the external momenta. Finally, using the external 
momenta of Eq. (4.1), we computed the amplitude introducing a non-zero mass mj for the 
fermions in the loop 8 . The results are plotted in Fig. 3, for the three cases nif = 0.5 GeV, 
nif = 4.5 GeV and nif = 12 GeV. 

The code we prepared for producing the results presented in this section is written in 
FORTRAN 90. Even if we did not spend too much effort in optimizations, it can compute 
about 3 phase-space points per second, when working in double precision. All figures in 
this section are actually produced by using double precision, but, to perform a realistic 
integration, we still need quadruple precision, that slows down the speed by about a factor 
60. We are working in implementing the expansions presented in the previous section with 
the aim of being able to perform a stable integration over the full phase space, that is a 
"proof of concept" for any method. 

8 We used here the scalar one- loop functions provided by FF [27]. 
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Figure 3: Helicity configuration [+ H ] for the momenta of Eq. (4.1) for different values 

of the fermion mass in the loop: m/ = 0.5 GeV (diamond), nif — 4.5 GeV (gray box) and to/ = 12 
GeV (black dots). The continuous line is the result for the massless case. 

5. Conclusions 

Computing the massless QED amplitude for the reaction 2j — > 47, although still unob- 
served experimentally, is a very good exercise for checking new methods to calculate one- 
loop virtual corrections. Such a process posses all complications typical of any multi-leg 
final state, for example a non trivial tensorial structure, but also keeps, at the same time, 
enough simplicity such that compact analytical formulas can still be used as a benchmark. 
However, it is oversimplified in two respects. Firstly, the amplitude it is completely mass- 
less. Secondly, the amplitude is cut constructible, namely does not contain any rational 
part. 

In the most general case of one-loop calculations, the presence of both internal and 
external masses prevents from obtaining compact analytical expressions. Then one has to 
rely on other computational techniques. For example, it is known that cut-constructible 
amplitudes can be obtained through recursion relations. But, even then, the presence of 
rational parts usually requires a separate work. 

For such reasons, it would be highly advisable to have a method not restricted to 
massless theories, in which moreover both cut-constructible and rational parts can be 
treated at the same time. Such a method has been introduced recently in Ref. [17] and, in 
this paper, we applied it to the computation of the six-photon amplitude in QED, giving 
also results for the case with massive fermions in the loop. We also showed in detail how 
the rational part of any m-point one-loop amplitude is intimately connected with the form 
of the integrand of the amplitude. Once this integrand is numerically computable, cut- 
constructible and rational terms are easily obtained, at the same time, by solving the same 
system of linear equations. This is a peculiar property of our method, that we tested in the 
actual computation of the six-photon amplitude. In practice, we did not use the additional 
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information on its cut-constructibility and verified only a-posteriori that the intermediate 
rational parts, coming from all pieces separately, drop out in the final sum. 

Finally, we presented all relevant formulas needed to infer the rational parts from the 
integrand of any m-point loop functions, in the renormalizable gauges. 

In addition, we presented, by analyzing in detail the 2-point case, an idea to cure the 
numerical instabilities occurring at exceptional phase-space points, outlining a possible way 
to build up expansions around the zeroes of the Gram-determinants. 

Having been able to apply our method to the computation of the massive six-photon 
amplitude, we are confident that our method can be successfully used for a systematic and 
efficient computation of one-loop amplitudes relevant at LHC and ILC. 
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Appendices 

A. Computing the extra-integrals 

In this appendix, we compute the extra-integrals listed in Section 2. Since a contribution 
0(1) can only develop for non-negative dimensionality V, the integrand in the Feynman 
parameter integral is always polynomial. First we decompose the integration as follows 



then, after using Feynman parametrization and performing first the integral over d e \i and 
then that one over d 4 q, one derives, for the extra-integrals of Eqs. (2.14)-(2.18) 





(A.l) 




(A.2) 
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where 



/poo ° 
[da} s = da ---da s 5(1-^2 a j), X s = P s 2 + M 2 S , 
Jo 1=0 

s s 
3=0 3=0 

In the following, we compute, as an illustrative example, the first three integrals of Eq. (A. 2). 
The remaining two can be obtained analogously. We start by changing the integration vari- 
ables as follows: 

«i = P1P2 ■■■ Ps 

«2 = PlP2 ■ ■ ■ Ps-l(l - Ps) 
«3 = PlP2 ■ ■ ■ Ps-2(1 - Ps-l) 



a s = pi(l - p 2 ) 

a = (1 - pi) , (A.4) 



so that 



/ [da] s = dpi I dpi--- I dp s p { { 1] p { 2 s 2) • • • p s _i , (A.5) 
J Jo Jo Jo 



from which one trivially obtains the first integral 



4n;2( S -l)) = _^ 2 r(j_l) +0(e) _ 



(A.6) 



To compute the second integral an integration over (P s )^, in needed. Since the integrand 
is symmetric when interchanging all fcj , we concentrate on the coefficient of, say, ki . Since 

/ dpi dpi--- dp s p[ s_1) p 2 S_2) • ■ -Ps-iocihfj, 
Jo Jo Jo 

= k 1/x / dpi dpi--- dp s pi s) P2 S_1) • • • P 2 s-iPs 
Jo Jo Jo 

= ki u —^ — - , (A.7) 
the final result reads 

4S 2(S - 1)} = - 2 ^| 5>V + 0{e) . (A.8) 



' 3 



To compute the third integral we need to integrate over the product (P s )u(Ps)v- Once again, 
given the symmetry of the problem, we can focus on the two contributions proportional to 
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ki^ki,, and k\ fjL k2v The first one gives 

/ dpi I dpi--- I dp s p^ 1] p^~ 2) ■ ■ ■ p s -iajk lfl k lu 
Jo Jo Jo 

= fci„fci„ f dpi I dp 2 ■ ■ ■ I d Ps p { i +1) p { 2 ] ■ ■ ■ pI-ipI 
Jo Jo Jo 



— k\^k\ v 



r( s + 3) ' 



(A.9) 



and the second reads 



dpi dp 2 --- dp s p^~^ p 2 S_2) •■■ Ps-iotiaiki^v 
Jo Jo 

/ dpi I dp 2 --- dp s Pi +1) p 2 s) ■ ■ ■ P 3 s -ips(l ~ Ps) 
Jo Jo Jo 



kink 



1 



Summing up all of the possibilities one obtains 



(A.10) 



^»;2(.-i)) = -^vr 2 ^^ < 

r(s + 3) 



s ^ s s 



3=1 



3=1 



(A.ll) 



B. The general basis for the 2-point functions 



In this appendix, we solve the problem of reconstructing the coefficients of the 2-point part 
of the integrand of any amplitude 



A(q) 



N(q) 
D D 1 ' 



(B.l) 



by assuming N(q) at most quadratic in q and k\ = {jp\ — po) ^ 0. In particular also the 
case of vanishing k\ is included. First, we introduce a massless arbitrary 4-vector v, such 
that (v ■ k\) ^ 0, that we use to rewrite k\ in terms of two massless 4- vectors (we also take 
£ 2 = 0) 



k\ = 1 + av , 



giving 

7 = 2 (h ■ v) = 2 ( 

Then, we introduce two additional independent massless 4-vectors £7^ defined as 

ef} = <e\-f\v], e%=<v\-f\£\, 



k 2 

■ v) and a = — 

7 



(B.2) 



(B.3) 



(B.4) 
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for which one finds 



(£ 7 -e 8 ) = -2 7 , (B.5) 
and we decompose + p$ in the basis of k\, v, l 7 and £g 

<f = -pfi + yk!t + y v v» + &£(} + y a e%, (B.6) 
so that N(q) takes the form 

N(q) = b + 6 [(« + Po) • «] + ^00 [(g + Po) • + &11 [(« + Po) • 4] + 621 [(g + Po) • 4] 

+ Mfa+Po) -^] 2 + &22[(9 +P0) -4] 2 

+ M(g + po) -Mte + Po) ■«] 

+ M(9+Po)-4][(? + Po)-«] + 0(£>o) + 0(£)i). (B.7) 

Notice that, because of the identity 

2 ( g • fci) = D 1 - D + (di - d ) , with di = m? - pf , (B.8) 

any term proportional to [{q + po) ■ k\] either contributes to the constant term b or it is 
included in the terms C(Do,i) we are neglecting 9 . The same happens for the combination 

[(<Z + Po)-M(9 + Po)-4]- 

To be able to determine all of the coefficients appearing in Eq. (B.7), disentangling 
completely the contributions O(l?o,i), we look for a q that fulfill the requirement 

Do = D 1 = 0. (B.9) 

For a q written as in Eq. (B.6) this implies the system 



V7V8 = F. 



u 



d 1 -d - 2yk\ 
Vv = , (B.10) 

7 



where 



F y = -±-(ml-y (d 1 - do) + y 2 kl) . (B.ll) 

It is convenient to introduce two classes of solutions. In the first class, that we call q^ k , 
we take y fixed and choose 2/7 = ±e l7T / k . In the second class, that we call q'^, we take y 
fixed but choose y§ = ±e m / fc . The coefficients b, bu, 621 , &12 an d 622 can be obtained by 
evaluating Eq. (B.7) at the values 

9oi> 9ra> I03 ' ( B - 12 ) 

or 

„/± „/± „/± 



9oi > 902 . 903 • ( B - 13 ) 



9 We suppose to determine them at a later stage of the calculation. 
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In the first case, the coefficients read 



b = b , 61 = -2 7 fe 21 , b 2 = 4 7 2 6 22 , 6_i = -27F0&11 , 6-2 = 4 7 2 F 2 6 12 , (B.14) 



with 



b±i = 
bo = 
b±2 = 



--[T-( qi )±iT-(q 2 )] , 
T+( qi )+T+(q 2 ) 



60) 



]_ _ e =F2Mr/3 



and where 



N(q+) ± N(qZ) 



(B.15) 



(B.16) 



In the second case, one obtains instead 

6 = 6, 6 / 1 = -2 7 6 11 , & 2 = 4 7 2 6 12 , ^-1 = -2 7 F 6 2 i 

with 



6'_ 



4 7 2 F 2 6 22 , (B.17) 



•±2 



T+(q[)+T+(q' 2 ) 



T+(q[)-T+(q> 2 ) _ e ±2i«/3 {T+{q , 3) _ ^ 



I — e =F2i7r/3 ' 



and where 



N(q'+)±N(q' 



OA; / 



(B.18) 



(B.19) 



The reason why we have chosen two sets of solutions is that, in some special kinematical 
configurations, Fo can vanish. Therefore, numerical stable solutions are obtained by taking 
621 and 622 from Eq. (B.14), and b\\ and 612 from Eq. (B.17), while b is well defined in 
both cases. 

The coefficients bo and 600 can be determined, in terms of additional solutions of the 
kind qj^ and , by defining the combinations 

S(q) = N(q) -6-611 [(q + Po) ■ i 7 ] ~ hi [(? + Po) " 4] 
- b u [(q + Po) -ir} 2 -b 2 2[(q+Po) • 4] 2 , 



S(q+)+S(q_ 



U(\) = 



as the two solutions of the system 



Al/ 



U(X) 
U(a) 



2 



4 

0-7 (T 2 7 2 



bo 

&00 



(B.20) 



(B.21) 
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The determinant of the matrix above is always different form zero, for non vanishing A and 
a, when a 7^ A, so that numerical inaccuracies never occur. 

Finally, the two last coefficients 601 and 602 are determined, in terms of q^ k and q'J~ k , 
as solutions of the system 



( Z{q + Xk ) \ 



( -\~f 2 F x e- in / k -AyV 7 ^ \ 
y -cr7 2 e i7r / fe -<r-i 2 F a e- i ' K / k j 



( lm \ 
\ b 2 ) 



(B.22) 



where 



Z(q) = S(q) - b [(q + p ) ■ v] - b 00 [(q + p ) ■ v} 2 . (B.23) 
Once again one verifies that when, for example, k = 3 the system never becomes singular. 
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